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Abstract 



Background: Stochastic biochemical reaction networks are commonly modelled by the chemical master 
equation, and can be simulated as first order linear differential equations through a finite state projection. Due 
to the very high state space dimension of these equations, numerical simulations are computationally expensive. 
This is a particular problem for analysis tasks requiring repeated simulations for different parameter values. Such 
tasks are computationally expensive to the point of infeasibility with the chemical master equation. 
Results: In this article, we apply parametric model order reduction techniques in order to construct accurate 
low-dimensional parametric models of the chemical master equation. These surrogate models can be used in 
various parametric analysis task such as identifiability analysis, parameter estimation, or sensitivity analysis. As 
biological examples, we consider two models for gene regulation networks, a bistable switch and a network 
displaying stochastic oscillations. 

Conclusions: The results show that the parametric model reduction yields efficient models of stochastic 
biochemical reaction networks, and that these models can be useful for systems biology applications involving 
parametric analysis problems such as parameter exploration, optimization, estimation or sensitivity analysis. 



Background 

The chemical master equation (CME) is the most basic mathematical description of stochastic 
biomolecular reaction networks 1 1 , 2 . The CME is a generally infinite-dimensional linear differential 
equation. It characterizes the temporal development of the probabilities that the network is in any of its 
possible configurations, where the different configurations are characterized by the molecular copy numbers 
of the network's chemical species. 

Due to its infinite dimension, the CME is usually not directly solvable, not even with numerical methods. 
A recent breakthrough in the numerical treatment of the CME was the establishment of the finite state 
projection (FSP) method by Munsky and Khammash 3 . They showed that it is possible to compute a 
good approximation to the real solution by projecting the CME to a suitable finite subdomain of the 
network's state space, and solving the resulting finite-dimensional linear differential equation on that 
domain. Nevertheless, the FSP approach still yields very high-dimensional models which are 
computationally expensive to simulate, even for small biochemical networks. The efficient simulation of the 
CME is an area of active research, and recently other simulation methods have been developed that can 
also be used for larger networks 4,5. 

Despite this progress, the direct simulation of the CME remains a computational bottleneck for common 
model analysis tasks in systems biology. It is especially problematic for tasks which require the repeated 
simulation of the model using different parameter values, for example identifiability analysis, parameter 
estimation, or model sensitivity analysis. Thereby, while a single or a few evaluations of a CME model 
with the FSP or other approaches may still be computationally feasible, the necessity of many repeated 
simulations will quickly render higher- level analysis tasks infeasible. 

Mathematical methods that approximate the behaviour of a high-dimensional original model through a 
low-dimensional reduced model are a common way to deal with complex models. Especially for linear 
differential equations, model order reduction is a well established field and several methods to compute 
reduced order models are available |6 . Note that the step of generating a reduced model is usually 
computationally more expensive than a single or even a few simulations of the original high-dimensional 
model. But the simulation of the resulting reduced models is frequently orders of magnitude faster than 
the solution of the original model. So, model reduction is worth the effort if many repeated simulations are 
to be expected. Unfortunately, for analysis tasks which require the repeated model simulation with 
different parameters, classical model reduction methods are not helpful. With these methods, the reduced 
model depends on specific parameter values in the original model, and the reduction needs to be redone for 
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different parameter values. Thus, for the mentioned analysis tasks, the model reduction process would have 
to be repeated for each new parameter value, and no gain in computational efficiency would typically be 



are not so suitable for parametric analysis tasks. 

Fortunately, model reduction methods where parameters from the original model are retained as adjustable 
parameters also in the reduced model are now being developed. These methods allow to compute a 
reduced model which uses the same parameters as the original model, and where the reduced model can 



The purpose of this paper is to introduce the application of these parametric model reduction methods to 
finite-state approximations of the chemical master equation, and to show possible usage scenarios of such 
an approach. The structure is as follows. In the following section, we introduce some background and 
notation concerning the modelling of chemical reaction networks and parametric model order reduction. 
We also show how the parametric model order reduction methods can in fact be applied to the CME. 
Afterwards, we apply the reduction technique on two reaction network models and corresponding 
parametric analysis tasks. 

Theoretical basics and methods 

We start with some preparatory background on the chemical master equation (CME) and parametric 
model order reduction. This serves in particular to fix the notation used throughout the remainder of the 
article. Then the application of parametric model order reduction to the CME is introduced. 

The chemical master equation 

The structure of a biochemical reaction network is characterized completely by the list of involved species, 
denoted as Xi, X2 . . . , X^, and the list of reactions, denoted as 



where m is the number of reactions in the network, and the factors aij G No and (fij G No are the 
stoichiometric coefficients of the reactant and product species, respectively 12 . The net change in the 
amount of species i occuring through reaction j is given by 



possible. While classical model reduction techniques have been applied to the CME in the past [7^, they 



directly be simulated with any choice of parameter values [8-11 . 



n 



n 




(1) 



Nij = (fij — (Tij. 



(2) 
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Reversible reactions can always be written in the form ([T]) by splitting the forward and reverse path into 
two separate irreversible reactions. 

For a stochastic network model, the variables of interest are the probabilities that the network is in any of 
the possible states which are characterized by the molecular copy numbers of the individual species 
Xi,X2 . . . We denote the molecular copy number of Xi by [Xi] G Nq. Then, the state variables of the 
stochastic model are given by the real numbers 

pit, xi,X2, ...,Xn) = Prob([Xi] = xi, [X2] = X2, . . . , [Xn] = at time t), (3) 

for G No, i = 1, . . . , n. As a short-hand notation for (|3|, we write p(t, x), with x G Nq • 

The transitions from one state to another are determined by chemical reactions according to ([T]). The 

changes in the molecule numbers are described by the stoichiometric reaction vectors 

Vj = {N^j N2j ••• Nnj)^ elT". (4) 
To avoid needlessly complicated cases, we assume Vj / Vk for j ^ k. 

The probabilities of the network being in any of the possible states x evolve over time, and their evolution 
is governed by the chemical master equation (CME) as derived by 1 1 . From a given molecular state x, one 
can compute the propensity that reaction j takes place according to the law of mass action as 

Mx,0)=0,f[('''), (5) 

where 6 = {6j)JLi is the vector of reaction rate constants, which are model parameters depending on the 
physical properties of the molecules involved in the reactions. The propensities are related to the 
probability that reaction j will occur in a short time interval of length dt when the system is in state x: 

Prob(reaction j occurs in [t, t -\- dt] \ [X] = x) = Uj{x, 0)dt + o{dt). (6) 

Taking the possible transitions and the corresponding reaction propensities together yields the chemical 
master equation (CME), a linear differential equation where the variables are the probabilities that the 
system is in each of the possible molecular states x: 

d ^ 

-^Pit^x) = ^{i^j{x - Vj,e)p{t,x- Vj) - Vj{x,e)p{t,x)), (7) 

for X G Nq. The CME ^ is subject to an initial condition p{to,x) = po{x) for x G Nq. 

Despite being linear, the CME is hard to solve numerically. This is due to the problem that the state space 

is for most systems infinite- dimensional, since all possible states x G Nq of the reaction network ([T]) must in 
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general be considered. Instead of directly solving the CME ([7|), a number of alternative approaches to 
study the stochastic dynamics of biochemical reaction networks have been suggested. The most common 
approach is to generate a simulated realization of the stochastic process described by the reaction network 
([T]), using for example the Gillespie algorithm 13 . In this approach, the probabilities p{t^x) for the 
possible system states are obtained from many simulated realizations. However, since this requires a large 
number of realizations, it is computationally expensive. 

As a more direct approach, Munsky and Khammash 3^ have proposed the finite state projection (FSP), 
where the CME is solved on a finite subset of the state space. Here, this subset is denoted by 1], and is 
defined as 

n = {x^'^ |z = l,...,4cN^, (8) 

where the x*^*^ are the system states for which the probabilities are computed in the projected model. The 
underlying assumption is that the probabilities for other states will be very low on the time scale of 
interest — otherwise the FSP may not yield good approximations to the solution of the CME. In particular 
we assume the time interval of interest to be given by [0, T] for final time T > 0. The probabilities for the 
states x^*^ in Q are written in the vector P{t) approximating t) at the finite number of states Q: 

Pit) = {Pm),=,_a - e [0, l]'^. (9) 

The equation to be solved with the FSP approximation is 

^P{t)=A{e)P{t) forte(0,T) 

dt (10) 

where A{6) G R^><^ is the matrix of state transition propensities, and Pq = (po(^^*^))^_]^ d '^^ ^ vector of 
initial probabilities for the states in Vt. The elements of the matrix A{0) are computed as 



r=l 



Aii{0) 



i'r{x^^\ 0) if x^^^ = x^^^ ^ r = 1, . . . , m 
otherwise. 

We will frequently omit the parameter dependence of the solution (and other parametric quantities). 



(11) 



Hence the solution as abbreviation of P(t,^), of (10) is an approximation to the solution p{t,x) of 

the orginal CME on the domain Q. Munsky and Khammash [s] have also derived an upper bound on the 
error between the solution P{t) computed via the FSP, and the solution of the original CME x) on ft. 
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Here, we consider in addition an output vector y eW defined by 

y{t) = CP{t), (12) 

with C G R^^^. Examples for relevant outputs are the probability that the molecular copy numbers are in 
a certain domain Q C which is achieved by the row vector output matrix C defined by Q = 1 if 
x*^*^ G Cl^ otherwise Ci = 0, with p = 1, or the expected molecular copy numbers, given by 

d 

ye{t) = Y,x^''>Pi{t), (13) 

i=l 

i.e. C = {x^^\ . . . , x^^^) with p = n. 

The basic motivation for the model reduction presented here is that we are interested in parametric 



analysis of the model, where the model (10) has to be solved many times with different values for the 
parameters 0. Due to the typical high dimensions of the matrix A{0)^ already a single simulation is 
computationally expensive, and analysis tasks requiring many repeated simulations are often 
computationally infeasible. Thus, the primary goal is to derive a reduced model which is rapidly solvable 
and provides an approximation y{t) to the output y{t)^ potentially without any consideration of the 
original state vector P{t). 

Order reduction of parametric models 

Model order reduction of parametric problems is a very active research field in systems theory, engineering 



and applied mathematics. We refer to [8j[T0 , 11 and references therein for more information on the topic. 
Here, we apply the reduction technique for parametric problems presented in adopted to our notation. 
It is based on two biorthogonal global projection matrices V, ly G W^^^ with r <^ d and W^V = Id, where 
r is the dimension of the reduced model. The matrix V is assumed to span a space that approximates the 
system state variation for all parameters and times. The construction of such matrices will be detailed in 
the next subsection. 

The gain of computational efficiency in repeated simulations comes from a separation of the simulation 
task into a computationally expensive "offline" phase and a computationally cheap "online" phase. In the 
offline phase, suitable projection matrices V and W are computed without fixing specific parameter values. 
With the projection matrices, a reduced model with the same free parameters as the original model is 
computed. In the online phase, the reduced model is simulated with the actually chosen parameter values, 
which is typically several orders of magnitude faster than the simulation of the original model. For analysis 



6 



tasks with repeated simulations, only the online phase has to be repeated for different choices of the 
parameter values, yielding an overall gain in computational efficiency. 



Decomposition in parametric and non-parametric part 

The reduction technique assumes a separable parameter dependence of the full system matrices and the 
initial condition. This means, we assume that there exist a suitable small constant Qa ^ N, parameter 
independent components A^^^ G R^><^ and parameter dependent scalar coefficient functions ^^^i^) 
= 1, . . . , Qa such that 

Qa 

A{e) = ^^^l\e)A^^^ (14) 

q=l 

and similarly for the system matrix C and initial condition Pq. We assume that C V stems from some 
domain V C of admissible parameters. In the next step, the reduced component matrices and initial 
conditions are determined by 

4^] := W^A^^W, C]^] := Ct^V, P^'J := P^'K (15) 

for g = 1, . . . , Qa' The resulting quantities A^r\ cl?\ and Pq^J are r-dimensional vectors or matrices and 
independent of the high dimension d. The basis computation and the computation of these reduced system 
components is performed once and parameter-independently in the offline-phase. Then, in the online-phase, 
for any new parameter the reduced system matrices and the initial condition are assembled by 

Qa 

and similarly for Pro(^) and Cr{0). The low dimensional reduced system that remains to be solved is 

^Pr{t) = Ar{0)Pr{t) for t G (0, T) 

Pr{0)=Pro{0) (17) 
y{t) = Cr{0)Pr{t). 

From the reduced state Pr{t)-, an approximate state for the full system can be reconstructed at any desired 
time by P{t) = VPrit). Also the difference between the approximated output y(t) and the output y{t) of 
the original model can be bounded by so called error estimators. A-posteriori error bounds for the reduced 
systems as considered here are given in 
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Basis generation 

Different methods for the computation of the projection bases V and W exist. In systems theory, methods 
hke balanced truncation, Hankel-norm approximation or moment matching are apphed, that approximate 
the input-output behaviour of a hnear time-invariant system [g). The resulting reduced models can be 
applied for varying input signals. Extensions to parametric problems exist, e.g. [sjjTT]. As we do not have 
varying inputs in the problem studied here, we consider snapshot-based approaches to be more suitable. 
This means, the projection bases are constructed by solution snapshots, i.e. special solutions computed for 
selected parameter values. 

The generation of projection matrices V and W must be done in such a way, that they are globally well 
approximating the system states over the parameter and time domain. A possible way to achieve this is 
the POD-Greedy algorithm, which has been introduced in 14 and is meanwhile a standard procedure in 
reduced basis methods p!5{p!6] . The algorithm makes use of a repeated proper orthogonal decomposition 
(POD) of trajectories P : [0,T] R^, which for our purposes can be defined as 



POD{P) := arg min / \\P{t) - {v^ P{t))v\fdt. 

veRd,\\v\\ = lJo 



(18) 



Intuitively, POD{P) G M'^ is a state space vector representing the single dominant mode that minimizes 
the squared mean projection error. Computationally, this minimization task is solved by a reformulation as 
a suitable eigenvalue problem. Consider the correlation matrix C = P{t)P{t)^dt. Then, 
^* = POD{P) e is an eigenvector corresponding to the largest eigenvalue Xmax of C, i.e.. 



Cv* = Arnax^*- For additional theoretical and computational details on POD we refer to 17 18 . We 
further require a finite subset of parameters V train C that are used in the basis generation process. As 
error indicator A(6>, V) we use the projection error of the full system trajectory on the reduced space 
spanned by the orthonormal columns of y, i.e. 



A(^,y):= / \\P{t,0)-VV^P{t, 
Jo 



0)\fdt (19) 



The POD-Greedy procedure which is given in the pseudo-code below, starts with an arbitrary orthonormal 
initial basis Vnq G ]R^><^o q^j^^ performs an incremental basis extension. The algorithm repeatedly identifies 
the currently worst resolved parameter (a), orthogonalizes the corresponding full trajectory with the 
current reduced space (b), computes a POD of the error trajectory (c), and inserts the dominant mode into 
the basis (d). 

function V = F OD-Gieedy {Vtr ain, Vnq, Stoi) 



1. N := No 

2. while Sn -= max^^-Ptra^r. ^(^^ ^n) > £toi 

(a) r := argmax^e7>t_. A(l9,FAr) 

(b) :=P(t,^)-F^(Vi^P(t,r)) 

(c) ^^+1 := POI)(^) 

(d) Vat+i := [Vat, ^AT+i] 

(e) N '=N^1 

3. end while 



Note that the algorithm is implemented such that the simulation of the full model, yielding P(t^6) in (19), 
is only performed once for each in the training set Vtrain- 

For concluding the basis generation, we set W := V . This satisfies the biorthogonality condition 



W^V = Id^ as V has orthonormal columns by construction. In practice the time-integrals in (18) are 
realized by a finite sampling of the time interval. 

A theoretical underpinning for the POD-Greedy algorithm has recently been provided by the analysis of 



convergence rates 19 . This is based on the approximation-theoretical notion of the Kolmogorov n-width 
dni^) of a given set C M^, which quantifies how well the set can be approximated by arbitrary 
A/'-dimensional linear subspaces of R^. The convergence statement for the case of exponential convergence 
then can be summarized as follows: If the set of solutions := {P{t^6)\t G [0,T],^ G V} C is compact 
and has an exponentially decaying Kolmogorov n-width dN{J~') < Me~^^'^ for some M, a, a > and all 
G N, then the error sequence {sN)NeN generated by the POD-Greedy procedure (cf. the definition in 
Step 2. in the pseudo code) also decays with an exponential rate, €n ^ CMe~^^^ with suitable constants 
/3, c, C > depending on M, a, a. Thus, if the set of solutions can be approximated by linear subspaces with 
an exponentially decaying error term, then the POD-Greedy algorithm will in fact find an approximation 
with an exponentially decaying error term, though possibly with suboptimal parameters in the error bound. 
Extensions of the POD-Greedy algorithm exist, e.g. allowing more than one mode per extension step, 
performing adaptive parameter and time-interval partitioning, or enabling training-set 



adaptation 16 20 . 
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Reduced models of the parametrized chemical master equation 

In this section, we describe how to apply the reduction method for parametrized models presented in the 
previous section to FSP models for the chemical master equation. 

As discussed in the previous section, the first step in the proposed reduction method is a decomposition of 
the d-dimensional system matrix A{0) as in ([l4|. Such a decomposition is possible for the case of mass 
action reaction propensities, as defined in ([5|, or generalized mass action, as recently suggested for the 
chemical master equation 21 . In this case, the length of the parameter vector is equal to the number of 
reactions m, and we decompose A{0) into m terms as 

A{0) = 6>i + • • • + . (20) 

Hence, concerning the notation given before, we have Qa = components A^^'^ and coefficient functions 
'^^a(^) ~ ^Q' Each matrix A^^^ in this decomposition comes from just the transition propensities 
corresponding to reaction and is defined by 

n 

4' = -11(4'^)'^" 



k=l 



k=l 

otherwise. 



More generally, such a decomposition is also possible if reaction rate propensities can be decomposed into 
the product of two terms, with the first term depending on parameters only, and the second term on 
molecule numbers only. This case is for example encountered when the temperature-dependance of the 
reaction rate constant is relevant, and the temperature T is a variable parameter in the Arrhenius equation 
= Ae RT , Since the output matrix C and the initial condition Pq are usually not depending on 
parameters in this framework, a decomposition of C and Pq is not considered. 
The situation is more difficult for reaction propensities involving for example rational terms with 
parameters in the denominator. The denominator parameters can not be included in the reduced order 



model by the decomposition outlined in (20) and (21). If variations in these parameters are however not 



relevant to the planned analysis, then they can be set to their nominal value, and the decomposition can 
directly be done as described above. Alternatively, approximation steps can be performed, such as Taylor 
series expansion or empirical interpolation 22 , that generate an approximating parameter-separable 
expansion. 
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Table 1: List of reactions and reaction propensity functions for the follicle switch model 23 . Nominal 

parameter values are ki = 4, = 75, Mi 25, ui 0.01^, V2 75, M2 = 25, U2 0.01^. 

Reaction Stoichiometry Vj Propensity 

Production of Xi (1,0)^ ui{ki + j^^J 

Degradation of Xi ( — 1,0)^ UiXi 

Production of X2 (0, 1)^ Mj^^^J 

Degradation of X2 (0, —1)^ U2X2 

Results for exemplary applications in genetic switching and oscillations 

In this section, we present the study of two example networks with the proposed model reduction method. 
With these examples, the applicability of the reduced modeling approach especially for analysis tasks 
requiring repeated simulations with different parameter values is illustrated. The first network is a bistable 
genetic toggle switch, where cells may switch randomly between two states, based on the model in |23 . For 
this network, the problem of parameter estimation with a reduced model is studied. The second network is 
a second-order genetic oscillator, based on 24 , where we perform a sensitivity analysis over a wide 
parameter range. 



Parameter estimation in a genetic toggle switch model 

Network description 



The genetic toggle switch considered here is an ovarian follicle switch model from 23 . It is a system of two 
genes which activate each other. The switch is modelled as a reaction network with two species Xi, X2, 
representing the gene products. The network reactions are specified in Table [l] 

In [23], this network was shown to describe a bistable switch with two probability peaks, one close to 
^(o//) ^ (0,0)^ and the other close to x^^^) = (^1,^2)^. 

In the study 23 , only the lower probability peak was of interest. Here, we are interested in the transition 
of the system from x^^^^"^ to x^^'^\ Therefore, the system is truncated to a rectangle 
(1 := {0, . . . , 150} X {0, . . . , 150} such that x^'^^\x^'^^ G Ct^ yielding a good approximation in the finite 
state projection to the infinite-dimensional chemical master equation. 

The next step is to apply the decomposition of the matrix A{0) as described in the methods section. Note 
that A{0) for the switch network contains rational terms with the parameters Mi and M2. Considering 
these two parameters as fixed quantities, the truncated CME for the follicle switch can be written as 

P{t) = (i^i^W + + uiA^^^ + ^2^^^^ + U2A^^^)P{t), (22) 
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where A^^, i = 1, . . . , 5 are of dimension 151^ x 151^ = 22801 x 22801. 

As initial condition we choose a probabihty distributed over some lower states 

f for xi + X2 < 20 

KO,x) = <^ 210 (23) 

t otherwise. 

For the parametric model reduction, we consider only variations in the parameters ui and 1^2- These 
influence both the steady state level of gene activity in the on-state as well as the switching kinetics and 
are thus of high biological significance in the model. Hence we set e := {ui, U2Y ^ [0.005, 0.02]^ as the 
parametric domain V. As final time we choose T = 10^ which corresponds to a time range of 
approximately 19 years, i.e. about three times the half-life time of the off-state estimated in |23 . 
Some state plots from the simulation of the full model are shown in Figure [l] These snapshots clearly show 
the transition of the switch from the off-state with low values for xi and X2 to the on state with high 
values. The parameter influence is mainly reflected in the speed of the transition: for the parameter vector 
(iii, U2) = (0.005, 0.02) in the lower row, most of the probability is already arranged around the on-state 
at the end of the simulation time. In contrast, for the parameter vector {ui^ U2) = (0.05, 0.005) in the 
upper row, a significant portion of the probability is still located around the off- state at this time point. 
Also, the transition paths are different: in the first case, the values for X2 are lower than the values for xi 
during the transition, while in the second case, this relation is reversed. 

As typical simulation time for a single trajectory of the full system, we obtain 98.2 seconds on a IBM 
Lenovo 2.53 GHz Dual Core Laptop. 

Basis generation 

We generated a reduced basis with the POD-Greedy algorithm, where the training set was chosen as the 
vertices of a mesh with 9^ logarithmically equidistant parameter values over the parameter domain V. We 
set etoi = 10~^^ as target accuracy. We use the projection error as error measure, hence precompute the 81 
trajectories for construction of the reduced basis. As initial basis we set Nq = 1 and Vn^ := Pq using the 
parameter independent initial condition. 

The POD-Greedy algorithm produces a basis of 33 vectors and the overall computation of the reduced 
basis takes 7.9 hours, the dominating computation time being spent in the error evaluations and POD 
computations. Some of the resulting orthonormal basis vectors are illustrated in Figure [2j The error decay 
curve and the selected parameters in the parameter domain are illustrated in Figure [3j We nicely observe 
an exponential error decay, which indicates a parametric smoothness of the solution manifold, cf. the 
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Figure 1: Illustration of some solution snapshots P{t) of the switch model (22) for parameter vector (i^i, 1^2) = 
(0.05,0.005) (upper row) and (2x1,^x2) = (0.005,0.02) (lower row) at times t = 0, 2 • 10^ 5 • 10^, and 1 • 10^ 
from left to right. 
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Figure 2: Illustration of the first eight basis vectors for the switch model generated by the POD-Greedy 
algorithm. 

convergence rate statement given before for the POD-Greedy algorithm. The selected parameters seem to 
be located at the boundary of the parameter domain, indicating that the model behaviour in between can 
well be interpolated from the model behaviours along the boundary of the parameter domain. 
The final reduced model of dimension 33 can then be simulated in 0.135 seconds, corresponding to a 
computational speedup factor of more than 700. 

Parameter estimation 

We exemplify a possible application of the reduced order model in parameter estimation, where we assume 
that a distorted output y{t) as the expected values E[xi] is available from population- averaged 
measurements. The task is to estimate the parameter values ui and U2 from such a noisy measurement. 
The reference parameter is 6ref = {ui^U2) = (0.01,0.01)^, and, for the purpose of this example, the 
measured output is produced by simulating the original model with the reference parameter values and 
adding 5% relative random white noise n{t) sampled from a standard normal distribution, 
Vmeasif) '= y{t^Oref){^ + 0.05n(t)). An illustration of the reference output y(t^Oref) and the noisy signal 

Umeas 

(t) is given in the left of Figure [I] 
We want to recover the values of the parameters ui and U2 based on fitting the reduced parametric model's 
output y{t,6) to the measured output ymeas{t)- As is commonly done in parameter estimation, we 
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Figure 3: Illustration of the error decay during the POD-Greedy algorithm (right) applied to the switch 
model and the selected parameters (left) being a small subset of the 81 training parameter points. 




Figure 4: Parametric analysis results for the reduced switch model: Application of parametric reduced models 
for parametric analysis: Illustration of the clean and noisy signals y(t^Oref) and ymeas{t)^ respectively (left), 
the optimization target J{0) over the parameter domain (middle), interactive parameter exploration by a 
graphical user interface (right). 
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formulate a least squares cost function as 

J{0) = r{ymeas{t) - y{t, e)fdt, (24) 
^0 

and estimate the parameters by 

6 est = argmin J(6>). (25) 

eev 

In such an optimization problem, typically many forward simulations are required for adjusting y to the 
measurement. This is a particular beneficial scenario for reduced order models, as these simulations can be 
computed rapidly. 

In order to gain a deeper insight into the optimization problem ([25|, we plot the values of the error 
functional J{0) over the parameter domain (middle of Figure [4|. Using the reduced model, the 
computation of the required 21^ = 441 trajectories is realized in less than one minute. This would be a 
significant computational effort when using a non-reduced model. 

From the cost function plot, we observe a narrow area of parameters which seem to produce a similar 
output as the reference parameter Oref- This shows that the two model parameters are not simultaneously 
identifiable from the considered output, and indicates that there may exist a functional dependence 
between the parameters ui and U2 such that the model yields similar outputs y{t). 

Assuming a functional dependence of Ui and U2 we now consider the 1-dimensional optimization problem 
along the line U2 = ^2, re/ = 0.01. We would like to recover ui from the optimization problem. The 
corresponding value of the cost function is J (Oref) = 3330.68, indicating a significant contribution of the 
noise. This restricted optimization problem is well conditioned and the optimization with a standard active 
set algorithm by MATLAB's command fmincon yields the estimated parameter Oest •= ('^i,es^7 0.01) with 
ui^est = 0.0100204, using 27 evaluations of the cost function. This accounts to a relative error in the ui 
value of 0.204%, hence excellent recovery. We refrain from plotting the recovered output y{t,6est) as it is 
visually indiscriminable from the output in the left of Fig. |4|. Interestingly, the optimization target value 
J{^est) = 3329.56 implies JiOest) < J {Oref) ^ which may stem from a slight approximation error in the 
reduced model or from the effects of the measurement noise. 

The right plot in Fig. |4] illustrates another application of reduced parametric models: We incorporated the 
model in an interactive graphical user interface in RBmatlab^ a matlab package for model order reduction, 
available for download at www.inorepas.org. This allows interactive parameter variations and 
instantaneous simulation response. 
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Table 2: List of reactions and reaction propensity functions for the oscillator model adopted from 24]. 
Nominal parameter values are ki = 15 7, k2 = 0.2, ks = I7, /C4 = 10^, k^ = 100 7, kQ = 6.5, kj = 100 7, 

s s s s s 

s = 10. 

Reaction Stoichiometry vj Propensity 



kis^ 

Degradation of Xi (—1,0)^ ksxi 



Production of Xi (1,0)^ , ^ 



Production of X2 (0, 1)^ k^s + 

Degradation of X2 (0,-1)^ 

Sensitivity analysis in a stochastic oscillator 

Network description 

The second case study is built on a genetic oscillator model showing stochastic resonance, which was 
presented in 24 . The oscillator is based on a negative feedback loop between two genes with one gene 
having positive autoregulation. The oscillator is modelled as a reaction network with two species Xi, X2, 
representing the gene products. The network reactions are specified in Table [2] In the original model 
in [24], the dynamics were described as stochastic differential equation for the amounts of Xi and X2, 
coming from a Langevin approximation to the stochastic dynamics 12 . For the framework used in this 
paper, the dynamics have to be described directly by the underlying CME. To achieve this, we introduce 
the parameter s which maps the dimensionless state variables from 24 to actual molecule numbers as 
required for the CME. Thus, s is also a measure for the network's noise level: the higher s, the larger the 
molecule number that is considered, and the smaller the noise level will be. 

The network model in Table |2] shows oscillations only in a stochastic description. The deterministic model 
has a unique asymptotically stable equilibrium point, but in a stochastic model, fiuctuations may push the 
molecular numbers beyond a certain threshold, inducing a dynamical response along a slow manifold. 



which corresponds to one oscillatory period 24 . Depending on the noise level, such responses will be 
initiated more or less often, corresponding to a more or less regular oscillatory pattern. 
The system is truncated to the rectangle 1] := {0, . . . , 300} x {0, . . . , 300}, which contains the relevant 
system states for the parameter ranges of interest. 

Similarly as in the switch example, the reaction propensity expressions contain rational terms in the 
parameters 5, /c2, and kQ. These three cannot be decomposed directly, so we do the decomposition 
described in the methods section for the other five parameters only. With this decomposition, the 
truncated CME for the genetic oscillator can be written as 

P{t) = {kiA^^^ + ksA^^^ + k^A^^^ + k^A^""^ + k7A^^^)P{t), (26) 
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Figure 5: Illustration of some solution snapshots P{t) of the oscillator CME model for parameter values 
/c4 = 15 (upper row) and = 30 (lower row) at times t = 0, 0.2, 0.6, 6.0 from left to right. 



where A^^, i = 1, . . . , 5 are of dimension 301^ x 301^ = 90601 x 90601. The initial condition for ([26| is 
chosen as a uniform distribution over the rectangle {0, . . . , 50} x {0, . . . , 50}: 



, ^ for xi < 50, X2 < 50 
1 ^ 



otherwise. 



(27) 



The time scale of interest for the model in (26) is for < t < T = 6. At the end of the interval, the 



probability distribution seems to approach a steady state. 

Some state plots are given in Figure [5] One observes a significant effect of the parameter on the 
amplitude of the oscillations. The simulation time for the detailed model was in average 7.3 minutes on a 
Dell desktop computer with 3.2 GHz dual-core Intel 4 processor and 1 GB RAM, without including the 
computation time for the construction of the state transition matrix A{0). 



Basis generation 

For the basis generation, the parameter k^ was assumed to vary within the interval [10, 100]. A reduced 
basis with the POD-Greedy algorithm was computed from a training set of 30 logarithmically equidistant 
parameters over the parameter domain (Figure [8|. As in the switch example, the target accuracy was 
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Figure 6: First 20 basis vectors for the oscillator model 



chosen as Stoi = 10~^^, and the initial basis was chosen from the initial condition Vi := Pq- 

The POD-Greedy algorithm produces a basis of 109 vectors, with an overall computation time of 16.5 

hours on the hardware as in the previous subsection. The first 20 basis vectors are shown in Figure [6j It is 

apparent that several of the basis vectors are directly included in order to reproduce the different 

amplitudes of oscillations that will occur under variations of the parameter k^. The error decay curve is 

shown in Figure [7| displaying an exponential error decay as also observed for the switch example. 

With the reduced basis V G ]r90601x109^ construct a reduced parametric model for the CME of the 

oscillator as 

Pr{t) = {hA^^^ ^ Al'^)Pr{t) 

(28) 

Pr{0) = F^P(O), 

with Ai'] = FT^[3]y ^ ^109X109 ^[o] ^ V^{hA^^^ + ksA^^^ + k^A^^^ + krA^^^)V e Rl09xl09^ Nq^^ ^^^^ 

since only k^ has been varied in the reduction process, the other parameters are no longer present as 
parameters in the reduced model, but just take their nominal values. While the same basis V could be 
used to construct another reduced model where all parameters are retained, it is unlikely that this other 
model will be a good approximation of the original one for varying values of the other parameters. 
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Figure 7: Results of the POD-Greedy algorithm for the oscillator model: Error decay curve for the oscillator 
model 

Sensitivity analysis of the oscillation amplitude 

As an application of the reduced order parametric model obtained in the previous section, we study the 
variations of oscillatory amplitude over a parameter range. Specifically, we consider 200 equally spaced 
values for the parameter in the interval [12, 40] and compute the probability that the amount of X2 is 
larger than 100: 

Prob(x2 > 100) = ^(^'^)' (29) 

x:x2>100 

with T = 6 the final time of the simulation. The results are shown in Figure [8] and show a clear decay of 
oscillatory amplitude for increasing values of /C4. Due to the significant time savings from the reduced 
model, this sensitivity curve can be computed with a high resolution. 



To evaluate the quality of the reduced model, we also computed the probability (29) using the original 



model (26) at two points within the considered interval for the parameter k^. As shown in Figure [sj the 
results from the original model are in perfect agreement with the predictions from the reduced model at 
these points. Since the points at which the original model was evaluated in this experiment were not part 
of the training set (shown as triangles on the parameter axis in Figure [8|, this shows that it is in fact 
possible to extrapolate the reduced model to parameter values that were not used to construct the basis. 



Conclusions 

In this paper, we have introduced the application of parametric model reduction methods to finite-state 
approximations of the chemical master equation. We have also presented two case studies where these 
methods are applied to CME models of different networks in order to make parametric analysis tasks 
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Figure 8: Sensitivity analysis of oscillation amplitude over a parameter interval. Blue line shows oscillatory 
amplitude over the parameter predicted from the reduced model. Red dots are validation results from 
a simulation of the original model. Triangles on the parameter axis indicate parameter values which were 
used in the construction of the reduced basis. 



computationally efficient. By this, it has become clear that parametric model reduction methods are a very 
useful tool for the analysis of stochastic biochemical reaction network described by the CME. 
Especially analysis tasks where many repeated simulations of a network with different parameter values are 
required can profit significantly from parametric model reduction. This includes for example sensitivity 
analysis or parameter optimization tasks such as identifiability analysis or estimation. Moreover, the 
significant speedup of the simulation for the reduced model allows an interactive exploration of the 
network's dynamics within the parameter space within a suitable graphical user interface. 
This contribution is just a first step in the application of parametric model reduction methods to the CME. 
One particularly important aspect that we have not discussed here is the computation of error estimates 
for certifying that the simulation output of the reduced model is within some tolerance of the 
corresponding simulation output of the original model. To maintain computational efficiency, the error 
estimation should be done without actually simulating the original model. Error estimation methods have 
been developed for parametric model reduction of generic models 9 , but tighter estimates could likely be 
obtained by taking into account the special structure of the CME models. Recent work for example refined 



the previous generic error bounds for stable models 25 
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